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TECHNICAL PUBLICATION 


NUMERICAL DETERMINATION OF CRITICAL CONDITIONS FOR THERMAL IGNITION 

1. INTRODUCTION 


1.1 Background 

Ignition or thennal explosion of a combustible substance occurs when exothermic reactions evolve 
heat so rapidly that it is impossible to preserve a stable balance between heat production and the loss of 
heat to the surroundings. This is an essential feature of many technological devices that employ auxiliary 
heat sources as a means of accelerating internal heat generation and engendering a runaway increase in 
temperature. On the other hand, there are numerous industrial processes involving the production and stor- 
age of reactive materials in which self-heating effects can culminate in spontaneous combustion or explo- 
sive effects and the primary concern is avoiding the occurrence of potentially hazardous circumstances . 1 

The archetypal example of self-heating is a porous pile of material in which heat is internally gen- 
erated by atmospheric oxidation. If the excess heat in the pile can be transported and dissipated to the sur- 
roundings fast enough, an equilibrium or steady state can be safely established. Under certain conditions, 
however, the dissipation mechanism cannot keep pace with the self-heating rate, and spontaneous ignition 
or explosion will occur. These critical conditions depend on the size and shape of the pile, the assembly 
temperature of the material, and the ambient temperature of the surrounding environment. 

Historically, assessment and control of self-heating hazards have been mainly conducted on an 
empirical basis, as established through long years of experience in the handling and processing of suscep- 
tible materials and products. The contemporary trend, however, is towards the development of reliable 
mathematical models and quantitative predictive procedures. In fact, the basic theoretical construct may 
be erected in a rather straightforward manner using conservation principles and well established descrip- 
tions of the underlying chemical and physical processes. The resulting mathematical formulation is com- 
monly referred to as the reaction-diffusion equation: 

pc„ = V (kVT ) + H(T) , (1) 

ot 

where t is the time, T is the temperature, p is the density, c p is the specific heat, k is the thennal conductiv- 
ity, and ll( T) is the rate of heat production per unit volume at temperature T. This nonlinear partial differ- 
ential equation is a localized expression of the conservation of energy and implies that the rate of change 
of thermal energy within a unit volume element is equal to the net conduction heat transfer through the 
bounding surface plus the volumetric heat generation rate. 



For most ignition problems of practical importance, it is possible to incorporate certain simplifying 
assumptions which make the mathematics more tractable. These include the following: 

• Negligible reactant consumption and reactant diffusion (i.e., “zero order” reaction). 

• Constant thermal conductivity. 

• Arrhenius temperature dependence for the exothermic reaction rate. 

Enforcement of these assumptions yields the following working form for the reaction-diffusion equation: 

a ?L = V 2 T + ^e- E °l RT , (2) 

dt k 

where a=pc p /k is the thermal diffusivity, Q is the heat of reaction per unit mass (i.e., ‘exothermicity’), A 
is the pre-exponential factor in the Arrhenius reaction rate tenn, E is the activation energy, and R is the 
universal gas constant. When applied over a bounded region Q, energy conservation principles yield a 
generalized boundary condition on the smooth surface <5X2 

- = h ( T «- T >) • O) 

dn 

where d/dn is the outward nonnal derivative on c)Q , h is the convective heat transfer coefficient, T a is the 
ambient temperature of the surrounding environment, and T s is the material surface temperature. 

Given the shape and size of a bounded region, appropriate boundary conditions, and values for the 
fundamental material properties, the basic objective is to mathematically exploit the reaction-diffusion 
equation and determine the critical parameters and conditions leading to the onset of ignition or thermal 
explosion. In particular, we are concerned with predictions for the critical ambient temperature, which 
defines an external environmental constraint for safe ‘storage,’ and the critical initial temperature, which 
defines an internal constraint for safe ‘assembly.’ 

1.2 Scope and Objective 

From a mathematical perspective, there are two fundamental strategies for attacking the reac- 
tion-diffusion equation and determining critical conditions for thennal ignition. These are the stationary 
(steady-state) model and the nonstationary (transient) model. 

In the stationary model, the time-dependent term in equation (2) is neglected and steady-state solu- 
tions are sought for which heat losses exactly balance heat production. This approach assumes unlimited 
reactants and implies that either a small steady-state excess temperature will become established in the 
body or the temperature will increase indefinitely. The principal attraction of the stationary modeling 
approach is a reduction of the problem to more amenable ordinary differential equation fonn, which has 
facilitated the development and refinement of standard mathematical methods capable of accounting for 
internal spatial temperature distributions and producing reliable estimates for the critical ambient tempera- 
ture. Because the stationary model cannot account for time evolution, however, it has only limited effec- 
tiveness in the prediction of critical initial conditions. It is well known, for instance, that many fires have 
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resulted from the assembly of reactive material at too high an initial temperature even though the storage 
conditions were sub-critical on the basis of steady state theory. See, for example, Bowes 1 (who refers to 
this as thermal explosion of the second kind), Rivers et al., 2 and Smedley and Wake. 3 

The nonstationary model, on the other hand, retains the complete time-dependent form of the reac- 
tion-diffusion equation and evolves the full temperature history of the self-heating body. The drawback 
of this approach is the general need to resort to numerical analysis and the fact that a full development 
history must be computed for each initial/boundary condition of potential interest. The distinct advantage, 
however, is that it can fully account for the initial assembly conditions and is therefore able to provide 
accurate estimates for the critical initial temperature. Weber et al., 4 for instance, recently conducted a lim- 
ited computational study of the nonstationary model which clearly demonstrated that the practical critical 
assembly temperature may, under certain circumstances, be 5-10% lower than the critical temperature 
obtained from steady-state theory. This has profound industrial implications in the assessment of fire 
hazards and the definition of fire safety standards, for which a clear distinction must be drawn between 
the classic ‘storage’ problem, where only steady-state temperatures are important, and the less recognized 
‘assembly’ problem, where the initial temperature threshold for self-ignition is of vital concern. 

Physical situations involving dynamical boundary conditions, where the ambient temperature or 
surface heat flux has a known time-dependent variation, represent an additional class of problems that 
require a nonstationary mathematical treatment for the accurate prediction of critical ignition conditions. 
Such ‘dynamic regimes of ignition’ might include the behavior of a combustible material under the action 
of an igniter with a predefined heat deposition rate or a reactive industrial stockpile exposed to diurnal 
variations in ambient temperature. Clearly, the standard stationary model results for ‘static regimes of 
ignition’ are inapplicable to this class of problem, and nonstationary approaches are the only viable course 
of action at the present time. 

The central objective of this research is the development and use of numerical techniques (1) to 
investigate the nonstationary solution sets of the full time-dependent reaction-diffusion equation subject 
to a general convective boundary condition and (2) to determine the critical threshold that distinguishes 
between initial conditions that evolve to a low-temperature steady-state and those that evolve to a high- 
temperature steady-state attractor (i.e., thermally ignite). This numerical methodology is then used as the 
basis for the following three-pronged research program: (1) conduct a broad ranging numerical study of 
the ‘assembly’ problem using a generalized one-parameter power law for the initial temperature profile; 
(2) investigate the relationship between the shape of the critical initial temperature distribution and the 
corresponding spatial moments of its energy content integral and attempt to forge a fundamental conjec- 
ture governing this relation; and (3) investigate the effect of dynamic boundary conditions on the classic 
‘storage’ problem and use the results of the nonstationary model to lay the groundwork for the develop- 
ment of an approximate solution methodology based on adaptation of the standard stationary model. 

1.3 Dimensionless Formulations of the Reaction-Diffusion Equation 

Evaluation and analysis of the reaction-diffusion equation is facilitated by the introduction of 
dimensionless parameters. From a historical perspective, it is important to note the traditional grouping 
of dimensionless variables suggested by Frank-Kamenetskii in the first comprehensive analytical treat- 
ment of the stationary model. 5 This classical formulation has been highly influential over the years, since 
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it facilitates certain simplifying approximations in the nonlinear Arrhenius rate tenn, making it amenable 
to analytical attack. In fact, virtually all subsequent theoretical developments have utilized this standard 
form as a common starting point. From a modem perspective, however, an alternative formulation based 
on the Bumell-GrahamEagle-Gray-Wake variables provides less circuitous contact with the primitive 
variables and enables a more direct physical interpretation. 6 Moreover, the facilitating features of the 
Frank-Kamenetskii variables are not necessarily advantageous for comprehensive numerical analyses 
where there is less need for simplifying assumptions. Thus, the mathematical treatment in this work will 
be constructed exclusively on the modern variable formulation. For the sake of completeness, however, 
and to facilitate translation and comparison between the two frames of reference, both formulations are 
briefly outlined in the following subsections. 

1.3.1 Frank-Kamenetskii Variables 


The Fra nk -Kamenetskii grouping of dimensionless variables was originally introduced as a means 
of simplifying the temperature dependence of the Arrhenius reaction rate term to permit the construction of 
exact analytical solutions. Because these variables appeared in the pioneering theoretical developments of 
the field, however, they became a de facto standard and tended to penneate all subsequent developments, 
despite certain drawbacks in clarity and interpretation. Thus, basic knowledge of the Frank-Kamenetskii 
formulation is essential as a frame of reference for understanding previous theoretical work. 

Development of this formulation begins with the definition of parameters for the dimensionless 
ambient temperature and the dimensionless reactant temperature, 


£ = 



and 


T -T 

6~ a 


eT n 


(4) 


Note that both parameters include the ambient temperature as a scaling factor and are therefore 
coupled. Thus, 0 relates the dimensionless temperature rise (T—T ) in the material to the characteristic 
temperature (eT a ). From the above definitions, one may readily establish the identity 


Eg = 1 , 0 

RT £ l+£d ’ 


(5) 


which leads to a dimensionless expression for the temperature dependence of the Arrhenius reaction rate 
tenn, e“ 1/e e 0/ (l+£0). Because £~ l =EJRT a is in the range of 10-100 for most materials, the nonlinearity 
exp[0/(l +£0)] is convex for small positive 6 and concave for 0>(£-£ 2 )/ 2. 


We next rescale and normalize the spatial and time coordinates by defining the following dimen- 
sionless variables: 


£ = y and r=-' T , (6) 

L aL 1 

where the characteristic length L represents the half-width of a symmetrical bounded region and s is the 
physical spatial coordinate referenced to the axis of symmetry. Substitution of the above dimensionless 
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variables into the reaction-diffusion equations (2) and (3) yields the Frank-Kamenetskii formulation in 
conventional compact form: 


^ = V;6> + <S// (l+e0) in Q 


(7) 


and 


dO 

% 


+ Bid = 0 


on dQ . 


( 8 ) 


Here, Bi=hL!k is the dimensionless Biot number and 8 is a dimensionless eigenvalue parameter given by 


8 = J*- l 2 PQ A e ~E a /RT a (9) 

RTa k 

It is convenient at the outset to develop a generalized formulation that is applicable to all three 
principal centrosymmetric solids, commonly referred to as the Class A geometries; i.e., the slab, infinite 
cylinder, and sphere. This is accomplished by expanding the Laplacian operator in cartesian, cylindrical, 
and spherical coordinates and observing that these can all be represented in the parameterized fonn: 



d~6 n dO 

if + ?^’ 


( 10 ) 


where n serves as a geometry selection parameter. That is, n = 0, 1, or 2 for the slab, infinite cylinder, and 
sphere, respectively. Hence, the final working fonn of the reaction-diffusion equation for class A shapes 
may be written as 


and 


— _ + } 1 — + 0 , n 

*• af 


in Q 


( 11 ) 


\-Bi0 = 0 on dQ (12) 

The Frank-Kamenetskii arrangement of the reaction-diffusion equation, in the form above, con- 
tains no additional simplifying assumptions and is appropriate for rigorous analyses of stationary and 
nonstationary models of thermal ignition. In fact, it has been the commonly used fonnulation in almost 
all previous mathematical and computational studies. 

Exact analytical attacks on the stationary model, however, require further simplifications to the 
nonlinear Arrhenius reaction rate term, and it is of passing historical interest to note the classic Frank- 
Kamenetskii approximation valid when e«l and 6 is not large. Under this assumption, equation (5) 
takes the simpler form 


— . (13 ) 
RT 
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One should note that the distinct advantage of the Frank-Kamenetskii variables — indeed the funda- 
mental reason for their original introduction — is the simplification provided by this approximation. Other 
simplified representations of the Arrhenius term, such as the quadratic approximation of Boddington et 
ah, 7 have been exploited for the development of exact analytical solutions in the Frank-Kamenetskii vari- 
ables. 


Therefore, these variables have been particularly useful in the construction of analytical solutions 
for the stationary model, where we are interested in the behavior of 0 under variations in the eigenvalue 
parameter <5. In carrying out the exact analysis, one finds that there are two solution branches for 0 when 
§<§ cr : the lower branch being stable and the upper branch being unstable. When 8>8 cr , however, no solu- 
tions exist. Thus, <5 may be considered as a bifurcation parameter in the Frank-Kamenetskii formulation 
such that 8 cr represents the first limit point of the corresponding bifurcation diagram. From a physical 
perspective, 8 cr is identified with the onset of ignition, since for 8>8 cr no steady-state solution exists and 
the temperature of the body will rise in time without bound. 


1.3.2 Burnell-GrahamEagle-Gray-Wake Variables 

Despite its predominance in the literature, the traditional Frank-Kamenetskii grouping of dimen- 
sionless variables has the effect of confusing the role of the ambient temperature when, in fact, it is the 
most practically significant control parameter in the problem definition. To circumvent this difficulty, 
Burnell et al. suggested an alternative dimensionless grouping with temperature rescaled independently of 
the ambient temperature. 6 

The major distinction in the new grouping is the definition of a dimensionless reactant temperature 
and a dimensionless ambient temperature that are completely decoupled: 


u = — and U = ^- . (14) 

E a E a 

As in the Frank-Kamenetskii fonnulation, we rescale and normalize the spatial and time coor- 
dinates using the previously defined dimensionless variables % = s/L and T =t/aL 2 . Then, substitution 
of the dimensionless parameters into the reaction-diffusion equations (2) and (3) yields the Burnell- 
GrahamEagle-Gray-Wake formulation in compact form: 


and 


— = Vju + ke~ 1 / u in Q 
dt ^ 


— + Bi(u-U) = 0 on 

88 , 


Here, A is a new dimensionless eigenvalue parameter given by 


(15) 


(16) 


6 



( 17 ) 


A= L 2 pQAR 
kE a 

In contrast with the eigenvalue S, A is decoupled from the ambient temperature. Thus, U is the only param- 
eter dependent on T a . 

The new formulation may also be put in a generalized form applicable to all three principal cen- 
trosymmetric solids. As before, we expand the Laplacian operator in cartesian, cylindrical, and spherical 
coordinates and observe that these can all be represented in the parameterized form: 



d~u n du 


(18) 


where n is the previously defined geometry selection parameter. That is, n = 0, 1, or 2 for the slab, infinite 
cylinder, and sphere, respectively. Hence, the final working form of the reaction-diffusion equation for 
class A shapes may be written as 


and 


du 

dr 


d~u n du . _i/,. 

+ Ae ' h 


% 


2 + 'z% 


in Q 


(19) 


1- Bi(u — U)=0 on dQ . (20) 

dz 

This arrangement of the reaction-diffusion equation contains no simplifying approximations 
beyond those previously contained in equations (2) and (3) and is a mathematically equivalent framework 
for the analysis of stationary and nonstationary models of thermal ignition. Indeed, for numerical analyses, 
there is no distinct advantage in using the Frank-Kamenetskii variables and it is indeed preferable to uti- 
lize the Bumell-GrahamEagle-Gray-Wake formulation, which provides a more direct link to the physical 
domain. Thus, all formal developments in this work are based on the latter formulation. 

In this new variable formulation, the most appropriate bifurcation control parameter is the dimen- 
sionless ambient temperature U with the dimensionless reactant temperature u serving as the response 
function. That is, A is normally fixed and the objective is to determine solution branches in the (u,U) plane. 
In this case, the first bifurcation point on the minimal branch of solutions occurs at the critical ambient 
temperature, U cr Thus, the critical ambient temperature can be inferred directly from the (u,U) bifurca- 
tion diagram as 


U 


RT, 


crit 


a,crit 


( 21 ) 


Table 1 summarizes the mathematical relationship between the (e, 9, 6 ) set of variables and the 
(u, U, 8 ) set of variables as an aid to interpretation and translation between the two formulations. 
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Table 1. Mathematical relationship between alternative dimensionless variables. 


(e, e, 8) (//, U, A) 

(//, U, X) (e, 0, 5) 

u = e(l+e6) 
U=e 

X= 8e 2 e l,E 

e=(u-u)/u 2 

e=U 

8= Xe~ l,u /U 2 


1.3.3 Shape Factor Method for Arbitrary Convex Regions 


The basic conceptual idea of using a geometric selection parameter n may be generalized to incor- 
porate arbitrary non class A three dimensional shapes by means of the shape factor method formalized by 
Boddington et ah, 8 which may be regarded as an extension of the equivalent sphere concept introduced by 
Wake and Walker. 9 For arbitrary convex body shapes of volume V and surface area S, the method produces 
non-integral values of n in tenns of the Semenov radius R s and the harmonic mean radius R 0 : 

r2 

,1 + 1 = 3 J-’ ( 22 ) 

where 


R 


s 



and 


( 23 ) 


R 


-ill 


dco 

~^2 

a 


( 24 ) 


Here, da is the solid angle subtended at the center O, and a is the radius from 0 to a point on the surface. 
The dimensionless spatial variable is defined in the conventional way, 


£ = 


a 

«0 


and, for class A geometries, L is related to n and R 0 by the fonnula: 


L - ciq = 



( 25 ) 


( 26 ) 
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Table 2. Values of R s , R {) , and n for some common geometries. 


Shape 

R s 

R o 

n 

Infinite slab 

3.000 

3.000 

0.000 

Rectangular parallelepiped (1:10:10) 

2.500 

1.731 

0.438 

Infinite cylinder 

1.500 

1.225 

1.000 

Infinite square rod 

1.500 

1.354 

1.443 

Rectangular parallelepiped (1:1:10) 

1.429 

1.354 

1.694 

Sphere 

1.000 

1.000 

2.000 

Equicylinder 

1.000 

1.115 

2.728 

Cube 

1.000 

1.194 

3.280 

Regular tetrahedral 

0.408 

0.537 

4.187 


A summary of results for class A and non class A shapes, as computed by Boddington et al, 8 is 
provided in table 2. This method has been extended to the revised variable formulation and applied to 
an extensive study of the stationary model by Balakrishnan. 10 This was accomplished by applying path- 
following techniques to the stationary form of equation (19) with the eigenvalue parameter X. The value 
of X for non class A geometries of unit size is recovered by scaling X, using the relation 


X = 


3 

(n + l)/?o 


X . 


(27) 


It is believed that this method is equally applicable to the nonstationary model and might prove 
practical for the estimation of critical initial conditions in assembly problems. 
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2. STATIONARY MODEL 


Historically, greater emphasis has been placed on theoretical understanding and mathematical 
modeling of the stationary model than the nonstationary model, because it reduces the determination of 
explosion or thennal ignition to a consideration of possible steady-state solutions of a standard nonlinear 
eigenvalue problem. These types of second-order ordinary differential equations have been extensively 
studied, both mathematically and computationally, and excellent reviews and compilations are available in 
the literature. 11 ’ 12 Such studies ordinarily adopt the F ra n k -K am en etski i variables and examine the behav- 
ior of 0 using 8 as a bifurcation parameter. Consequently, attention has been focused almost exclusively on 
detennination of 8 cr , the first bifurcation point on the minimal branch, which has direct practical implica- 
tions to the storage problem but only limited utility, at best, to the assembly problem. 

Balakrishnan, in one of the first attempts to utilize the Burnell-GrahamEagle-Gray-Wake variables, 
revisited the classic stationary model and carried out extensive numerical analyses of the eigenvalue prob- 
lem using path-following techniques. 10 This work illuminated the characteristic branch structure in the 
(u,U) plane, including accurate identification of U cr , and convincingly demonstrated that these alternative 
variables provide improved physical clarity and more straightforward interpretation in terms of the ambi- 
ent temperature control parameter. 

The major focus of this thesis, following the initial lead of Weber et al., 4 will be the solution of the 
nonstationary model in alternative variable form and determination of critical initial conditions in the oft 
neglected but practically important assembly problem. In doing so, it will first be necessary to construct 
stationary model solutions to provide a frame of reference for proper interpretation of the computed criti- 
cal initial condition bifurcation branches. Thus, a numerical analysis procedure is developed herein for the 
standard stationary problem in alternative variable form. Rather than follow the sophisticated path follow- 
ing method outlined by Balakrishnan, 10 however, a more direct and less laborious route is taken whereby 
the conventional two-point boundary value problem is reformulated as an equivalent initial value prob- 
lem. Before embarking on this development for the stationary model, it will prove useful to first review the 
qualitative solution behavior for the general eigenvalue problem as obtained from the time-independent 
form of equation (15): 


V\u + Xe 1 u =0 in Q 

where the boundary condition defined by equation (16) still applies. 


(28) 


2.1 Qualitative Structure of Solution Branches 

Consider the anticipated functional behavior of heat generation and heat loss in a porous pile of 
material undergoing chemical oxidation within a fixed temperature ambient environment. First, we note 
that the heat loss rate is assumed to be proportional to the temperature difference between the body and 
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the ambient environment and would therefore have a linear functional dependence on body temperature, 
as illustrated by the dashed lines in figure 1 corresponding to T a <T acr , T a = T acr and T a >T acr . On the 
other hand, we clearly expect the heat generation rate to vary in a nonlinear way with increasing body 
temperature. At relatively low body temperature when the driving potential for heat loss is low, the heat 
generation would tend to increase strongly with any incremental rise in temperature, whereas at high body 
temperature, when the driving potential for heat loss is higher and diffusion becomes a controlling factor, 
the heat generation rate would exhibit a weaker response to any incremental increase in temperature. Thus, 
the functional dependence of heat generation rate on body temperature would tend to follow a sigmoid 
type curve as illustrated by the solid curve in figure 1 . The intersection points of the heat generation and 
loss curves represent the steady state solutions we seek for the stationary storage problem. 



Figure 1 . Heat balance characteristics for the stationary ignition model. 


When T a < T u , the body is able to efficiently eject heat to the environment at a relatively high rate, 
and we obtain three steady-state solutions at intersection points 1,3, and 4. For this ambient condition, 
the stationary model implies that the body temperature will always approach 7j if 7j< Tj and T 4 if T [ > 7k . 
Here, T- is the initial assembly temperature. The intermediate intersection point at 7k is an unstable steady 
state and cannot be realized. When T a = T a cr the body ejects heat at a reduced rate, and we now obtain 
two steady-state solutions at intersection points 2 and 5. At this critical ambient temperature, the station- 
ary model implies that the body temperature will always approach T 2 whenever 7j< T 2 and T 5 whenever 
Tj> T 2 . When T a > T a cr , however, the ability to eject heat is significantly curtailed, and we obtain a single 
steady-state solution at intersection point 6. Thus, the body temperature will always approach T 6 for any 
initial value of 7], In most practical cases, the upper branch of steady-state solutions corresponds to tem- 
peratures in excess of the material flame temperature, and ignition will occur before those conditions can 
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be attained. Therefore, the upper branch serves as a nonphysical mathematical attractor only. In reality, 
the stationary model neglects all transient effects and is only applicable to storage problems where T i is 
well below the intennediate steady-state solution. It cannot be used to reliably predict a safe assembly 
temperature. 

There are certain mathematically significant conclusions of fundamental importance to be noted 
from this qualitative inspection of the stationary reaction-diffusion equation: 

• At least one steady-state solution exists for any finite value of physical parameters when T a > 0. 

• A multiplicity of steady-state solutions may exist whenever 0 <T a < T a cr . When T a < T a cr the body of 
reacting material will either tend to a stable lower steady-state temperature below the ignition point 
or to an upper steady state temperature above the ignition temperature. The critical material assembly 
temperature, 7} cr , demarking the two possible responses can only be roughly estimated using the unsta- 
ble intermediate steady state solution branch. Accurate detennination of T i cr requires consideration of 
the fully transient reaction-diffusion equation. 

• When T a = T a cr two steady state solutions exist. The lower solution corresponds to the first bifurca- 
tion point at the intersection point of the lower and intennediate steady-state branches, and the upper 
steady-state temperature exceeds the value required for ignition. 

• When T a > T a cr , only one solution exists and the body will always tend to the upper steady-state branch, 
in which case ignition is assured. 


The characteristic behavior of the steady-state solution is also dependent on variations in the value 
of A. 13 If A is small, a unique solution always exists. If A is sufficiently large, multiple solutions exist 
whenever T a < T a cr When T a > T a cr , a unique solution exists for any value of A. The transitional value of A 
above which multiplicity first occurs is denoted as A fr . When U=0, the transitional value of A above which 
multiplicity first occurs is denoted as A'. 

2.2 Numerical Methodology 

The possibilities for constructing analytical solutions to the stationary reaction-diffusion equation 
in the dimensionless set of variables are extremely limited for even the simplest shaped regions. Thus, 
numerical methods are required to obtain accurate solutions and reliable prediction of the critical param- 
eter values. The mathematical object of consideration is the solution set {( u,U , A)} to the second-order 
nonlinear ordinary differential equation for arbitrary shape factor n 

il + l*L +Xt -V. 0, 0<«<1 (29) 

df £d£ 

subject to the following boundary conditions at the axis of symmetry and the body surface: 
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du _q 
^ , E =n 


— +Bi(u s -U) = 0 , 
£=1 


where u s is the dimensionless temperature at the body surface. 

2.2.1 Reduction to First-Order ODE System 

The numerical solution of equation (29) is facilitated by reduction to a first-order system of dif- 
ferential equations. By defining u l =u and u 2 =du l /d^, for example, we obtain the first-order system 

du | 
dg 

duj n , -\/ u . 

— L- — u 2 -Xe V“i . 

Note that this system is nonautonomous since g appears explicitly in the denominator, and there is 
an apparent singularity at the origin, even though a solution must exist at that point. The singularity can be 
removed, however, by observing that when E, approaches zero, 

(32) 

£ £ £ 

since u{ (0) = 0. This approximation may then be used to eliminate u 2 in equation (3 1) as g -^0: 

= as (->0, (33) 

dg dg 


which gives 


du2 _ A _i/ Ml 
dg n + 1 


as E, — > 0 . 


Thus, we may redefine the first-order system, equation (3 1), in piecewise form as 


A _i 


dui n + 1 


di -%2-Ae- 1 ^ , 


where £ is a finite value very close to zero. Transfonnation of the boundary conditions defined by equation 
(30) yields the relations: 
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u 1 (l) = U-u 2 {l)/Bi 


and 


u 2 (O) = 0 . 


(36) 


2.2.2 Disjointedness and Redefinition 

The first-order system, equation (35), with boundary conditions, equation (36), defines a nonlin- 
ear two-point boundary value problem which may be solved using various sophisticated path-following 
methods. A serious numerical difficulty arises, however, when the value of A increases beyond A' and the 
extinction point — i.e., the upper limit point in the bifurcation diagram — is moved into the physically unac- 
ceptable domain of negative U values . 13 ’ 14 In this case, the upper and lower branches of solutions become 
disjoint in the physically acceptable region where U> 0. Furthermore, Balakrishnan has demonstrated that 
there can be no solutions with bounded derivatives in the case U< 0 for which u(r) changes sign in the 
region 0<£< 1. 10 

A straightforward approach for making a connection between the upper and lower branches of 
solution curves in the U< 0 region is to simply switch off the source term. 10 Mathematically, this may be 
formally accomplished through inclusion of the Heaviside function, H(u), as a multiplicative factor to the 
source term, in which case, 


where 


,2 

a u 


d<r <= 


(37) 


H(u) = 1 
H(u) = 0 


u > 0 , 
u < 0 . 


(38) 


Thus, the disjointedness may be removed and the system made autonomous by introducing the 
new variable u 3 = q and redefining the piecewise first-order system, equation (31), as follows: 


du\ 

dq 


= u 2 


du 2 

dE, 


du 3 

dq 


A _ 


V«i 


77 + 1 


0 


-77 — - Ae _1 /“i 
i/ 3 

u 2 


-77- 


77 3 


with the boundary conditions 

u l (l) = U-u 2 (l)/Bi , 


“3 

and 

0 

A 

s 



and 

Ui — 0 

(39) 

u 3 >g 

and 

0 

A 

a 

u 3 >g 

and 

0 

VI 

s' 


= 0 , 

and 

77 3 (°)=0 . 

(40) 
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2.2.3 Path-Following Method 

A numerical procedure is now needed to path follow the solution curves for the stationary thermal 
ignition model. Standard path-following techniques (viz., piecewise linear or predictor-corrector methods) 
are available to handle this problem, 15-22 but we shall develop a more straightforward approach based on 
direct numerical integration of the first-order system, equation (39), with boundary conditions, equation 
(40). Although this procedure will lack the sophistication of standard methodologies, it will be seen to 
be an efficient and highly effective means for constructing accurate bifurcation branch diagrams for the 
steady-state ignition problem. 

The central technical difficulty is that we are faced with a two-point boundary value problem, since 
1^(0) = i/(0), the dimensionless center-point temperature, is not known a priori. Otherwise, the system of 
equations could be readily solved as an initial value problem using conventional numerical integration 
techniques. Therefore, a simple ad hoc procedure has been adapted in which the solution branches in the 
{u(0),U} plane, with U as the bifurcation parameter, can be constructed from the solution of an equivalent 
initial value problem. 

The procedure is as follows: 

• Take u Q as an independent parameter for the dimensionless center-point body temperature such that 

0 — u O — u O,ma\ ■ (41) 

• Specify a value for A and use a standard fourth-order Runge-Kutta algorithm to numerically solve sys- 
tem equation (39) as an initial value problem from £=0 — » 1 with the initial conditions defined by 

ui(0) = u o , m 2 (0) = 0 ’ an d 1/3(0) = 0 . (42) 

• Use the computed boundary values /^(l) and i/ 2 (l) from the preceding step along with the specified 
Biot number to compute the corresponding body surface temperature from the boundary condition 

U = u l (l) + u 2 (l)/Bi . (43) 

• Repeat the procedure above for incremental values of u Q over the range 0 — > z/ 0 max and path follow the 
solution branch u(0) = u o versus the bifurcation parameter U. 

2.3 Validation 

An extensive numerical investigation of the stationary model was previously carried out and thor- 
oughly documented by Balakrishnan using conventional path-following tools, 10 and his results are now 
adopted as a benchmark validation standard for the ad hoc path-following procedure described above. 
Thus, the solution sets {(//,U,A)} are examined and validated for the principal centrosymmetric solids 
(n = 0, 1, 2), keeping U as the bifurcation parameter and allowing A and Bi to vary. 

For illustrative purposes, representative bifurcation diagrams were computed using the above path- 
following procedure for various values of A, assuming an infinite Biot number. These results are shown in 
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figures 2, 3, and 4 for n = 0, 1,2, respectively. It is found that these predictions are in excellent agreement 
with Balakrishnan’s calculations, which were obtained using more sophisticated methods and standard 
algorithms. 23 Both methods yield identical values of X lr and A' for all three geometric shapes, as indicated 
in the accompanying figures, and the corresponding critical dimensionless ambient temperatures U cr are 
in complete accord. 



u 

Figure 2. Computed bifurcation diagram for the infinite slab (n = 0) with Bi — > 


As a further validation step, comparative calculations were completed for U cr and w(0) for the prin- 
cipal centrosymmetric solids (n = 0, 1, 2), keeping U as the bifurcation parameter and allowing A and Bi 
to vary. The variation in the critical parameters against A are shown in figures 5 and 6 under the constraint 
Bi — > oo. The variation in the critical parameters against Bi are shown in figures 5 and 6 under the constraint 
A= 50. These results are all in excellent agreement with previously established theoretical predictions. For 
infinite Bi, it is observed that U cr and «(()) both decrease at the ignition point with increasing A. For a fixed 
value of A, we observe that U cr and //(()) both increase with increasing Bi and asymptotically approach 
limiting values as Bi — » °° (figures 7 and 8). 

Because the predictions were found to be precise and accurate, under all circumstances, in com- 
parison to established theoretical calculations, it is asserted that the results from the proposed ad hoc 
path-following method are fully valid and may be used as a reliable benchmark standard for the stationary 
model. Thus, results of the theoretical investigations to follow may be referred to this baseline with a high 
degree of confidence. 
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In most practical cases, the value of X is orders of magnitude higher than the typical range associ- 
ated with X tr and X'. Flammable materials of common interest, for instance, will normally fall in some 
range around X ~ 10 6 . Therefore, the bifurcation diagrams for the principal centrosymmetric solids have 
been constructed for X= 10 4 , X= 10 6 , and X= 10 8 with Bi -> and these results, summarized in figures 
9-11, will be used as a general frame of reference in the theoretical investigations of the nonstationary 
model to follow. 



Figure 9. Computed bifurcation diagrams with X = 10 4 and Bi — > 
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3. NONSTATIONARY MODEL 


Real world problems, such as the safe assemblage of hot reactive materials and so-called ‘hot 
spot’ ignition, where localized high temperature regions expand and initiate thermal explosion, can be 
significantly dependent on the contribution of transient heating processes. Thus, the fate of a self-heating 
material depends upon circumstances of assembly and external environment, and the most indiscriminate 
case requires consideration of the fully transient reaction-diffusion equation. 

Exact analytical constructions for the generalized nonstationary model have proved evasive, how- 
ever, and progress along these lines has been mainly limited to less general problem fonnulations derived 
from simplifications to the nonlinear Arrhenius rate term. The elementary treatment of Frank-Kamenetskii, 
for example, reveals salient features of the solution, 5 and Gray and Harper have developed more exact 
analytical representations using a quadratic approximation for the Arrhenius temperature dependence. 24 ’ 25 
More rigorous analytical work using expansion procedures in the region £<<1 has provided accurate 
self-heating and explosion solutions over a time span ranging from initiation to completion. 26-29 For a 
complete comprehensive treatment valid over the entire parameter range, however, numerical integration 
techniques are required. 

As recourse, we embrace a numerical approach to the Bumell-GrahamEagle-Gray-Wake form of 
the transient reaction-diffusion problem for class A shapes, as derived in section 1 and repeated here for 
convenience: 


du d"u n du 

dt d% 2 Z % 


in Q 


(44) 


— + Bi(u-U) = 0 on dQ , 

where we have introduced the parameter S= to represent the source term. 


(45) 


In proceeding, it is recognized that the reliability of computational solutions to partial differential 
equations is highly dependent on the integrity of the numerical scheme and the attention to detail and 
degree of care associated with its implementation. Therefore, the numerical methodology and implemen- 
tation procedures to be used in this study are thoroughly explained, developed, and validated as a major 
point of departure. 


3.1 Numerical Methodology 

A numerical solution of the transient reaction-diffusion equation, which is parabolic in time and 
elliptic in space, consists of a finite set of numbers from which the spatial distribution of the dependent 
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variable can be constructed at some instant in time. Thus, a numerical discretization methodology rests on 
the construction of a set of linear algebraic equations for unknown dependent variable values at a finite 
number of spatial locations — i.e., nodes or grid points — and the prescription of an algorithm for solving 
the algebraic system and advancing the solution over discrete increments of time. 

Systematic spatial discretization is accomplished by introducing finite difference expressions 
whereby the value of the dependent variable at each node is related to the values at a small number of 
neighboring grid points, only. The number of neighboring values included in the finite difference expres- 
sion, commonly referred to as the nodal support, detennines the numerical accuracy of the discrete spatial 
derivative. A common systematic approach for temporal discretization, which relies on the one-way trans- 
mission of infonnation in time, is to simply relate the initial dependent variable distribution to an evolved 
distribution over some small positive time increment. 


3.1.1 Generalized Discretization Equation 


3.1. 1.1 Interior Nodes. Since very high order accuracy is not essential for this problem, the spa- 
tial derivatives in the dimensionless reaction-diffusion equation are approximated as three-node support 
finite difference expressions via truncated Taylor series. For equally spaced grid points, the Taylor series 
expansion about any node i may be truncated after the third tenn to yield the following well known central 
difference approximations exhibiting second-order spatial accuracy: 


du 

x, 


u i + 1 u i - 1 

2AE, 



(46) 


d 2 u 


u i + 1 + u i - 1 

AE 2 



(47) 


Introduction of these relations into the continuum equation and temporal discretization over a 
small time increment At produces a generalized discretization equation centered on each interior node i: 
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_ m 
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2A<f; 


+fs, 

+(i -f)s, 


m 


(48) 


Note that subscripts and superscripts denote spatial and temporal indices, respectively, and that /is 
a weighting factor (0 </ < 1), which determines the relative influence of initial and final time step values 
during the temporal evolution process. 


For certain specific values of the weighting factor f the discretization equation reduces to one of 
the well known integration schemes for parabolic partial differential equations, as summarized in table 3. 
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Table 3. Relationship between weighting factor and integration scheme. 


Weighting Factor,/ 

Scheme 

0.0 

Fully explicit 

0.5 

C rank-N icholson 

1.0 

Fully implicit 


The fully explicit scheme expresses the value of the unknown dependent variable at node i and 
future time step m + 1 explicitly in terms of known neighboring values at the initial time step m . Conversely, 
the fully implicit scheme expresses the value of the unknown dependent variable at node i and future time 
step m + 1 implicitly in terms of the unknown neighboring values at the future time step m+ 1. As a middle 
path, the Crank-Nicholson scheme expresses the value of the unknown dependent variable at node i and 
future time step m + 1 partly in tenns of known neighboring values at the initial time step m and partly in 
tenns of the unknown neighboring values at the future time step m + 1 . 

It should be pointed out that the fully explicit scheme is prone to numerical instabilities and requires 
the utilization of extremely small time increments to obtain physically realistic results. It is, for all prac- 
tical purposes, not useful for serious calculations. The fully implicit scheme, on the other hand, has the 
distinct advantage of being unconditionally stable and will yield physically accurate results over relatively 
large time increments. Thus, it is often utilized as a means of incorporating computational robustness. For 
best accuracy, the Crank-Nicholson scheme is the superior choice, provided the time increment remains 
relatively small. 

Here, the discretization equations will be developed in the most general possible fonn with the 
weighting factor retained as an arbitrary parameter such that the time integration scheme may be selected 
as a matter of choice in the study. For practical use, it is advantageous to multiply the generalized discreti- 
zation equation (48) by At; and regroup common terms: 
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(50) 
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This result applies to all interior nodes (i=2,..., N- 1). 

3. 1.1.2. Boundary Nodes. Closure of the computational domain requires specification of bound- 
ary conditions at the symmetry axis and at the external surface c)Q and construction of discretization equa- 
tions valid at the corresponding boundary nodes. 

The development for the symmetry boundary condition, duld^ | ^_ () =0. is straightforward and 
directly follows from a forward difference approximation, 


du 

% 


< 3=0 


U2 - U\ 


+ o(A%) = 0 , 


(51) 


from which we deduce the desired discretization equation at the symmetry node (/= 1), 


m + 1 
U\ 


m + 1 
u 2 


(52) 


As a first step in the development of a discretization equation for the surface boundary node, con- 
sider a 3 -node support centered at i=N where a fictitious nodal value u* has been introduced as illustrated 
in figure 12. Then, introduction of the central difference approximation for the derivative in the boundary 
condition, equation (45), gives the expression: 


U N -i 

o- 




t/* 

o — ► 

Fictitious 

Node 


Figure 12. Computational grid structure at external boundary, including fictitious node. 


du 

% 


u*-u N _i 


= Bi(u — U] y) 


l£=l 

which may be solved for the fictitious dimensionless temperature u*: 

u* = 2A^Bi(U -u N ) + u N _i . 


(53) 


(54) 


Forming the central difference approximation for the second spatial derivative at the surface boundary 
yields the result 


a 2 

d u 
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u * +u j\f—i — 2 u ]\f 

A? 


( 55 ) 
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and elimination of u* using equation (54) provides the relationship 


d~u 2.A^BiU — 2(l + — 2n^y_j 

* <56) 

Now, we allow U to be a time dependent variable, introduce equations (53) and (56) into the reac- 
tion-diffusion equation (44), and discretize over a small time increment At to produce a generalized dis- 
cretization equation centered on the surface boundary node N: 

u 'l v +l ~ u N 2A%BiU m+1 -2(1 + A%Bi)ujy +l + 2u™l\ (nBi\t m+1 m+ \\ 

3 ? ~ UN > 

+ Ls m+1 +(i- f) 1^>BiU m -2(1 + A^Bi)u’^ + 2 u^_ x 
2 N y 11 A% 2 

Hi-f){^\{u m -u t S)+ i ^-S’S . (57) 

Note in this case that the source term has been halved, since the boundary cell volume external to 
the computational domain is half that of an interior cell. For practical purposes, we multiply equation (57) 
through by At, and regroup common terms to obtain: 


\ m+1 

m r 


U]y — f\ 2Bi + 


nBiAt, 


U m+l +(!-/) 2 Bi + 


nBiAt, 


, f ( ^ ^ m+1 , ft j:-\( ^ \ m e 2(l + A^Bi'j nBiAt, m+1 
+ f Te \ u N-l + \l-f) 77 \ u N-l~ f 77 + — 7 U N 


-(1 -/) 


2(l + A^Bi^ nBiAt, m c m+ 1 , t\^> c tn 

At + ~i7~ Un+/ Y Sn +,1_/) Y Sn 


A^ (\ + AE,Bi nA^Bi ^ m+1 _( 2 1 m+1 

= m UN -' 

, \ A % n r\( 1 + A ^ Bi , nA£Bi\ 1 m 

+ [TT- {l - f \-^r + ^r)\ UN 

+ 2 Bi + ^i- [fU m+l +(l-f)U m ] + f^S% +l +(l-f)^-S% . (59) 


26 



3.1. 1.3 Source Term Linearization. The objective of the discretization procedure is to reduce 
the partial differential equation problem to a set of linear algebraic equations, which may be effectively 
attacked using powerful linear algebra techniques. Thus, an additional linearization approximation must 
be introduced for the Arrhenius source term in the preceding development, since it contains a nonlinear 
exponential function, S=Xe l/u . 


The simplest approach in this case is to treat the source tenn in a fully explicit manner such that 


g>n+l _ gin 


i = l,...,N . 


(60) 


Utilization of this approximation is straightforward and requires no further explanation. 


An alternative and somewhat more sophisticated method relies on an expansion of the form 


(61) 


dS _ dS du 
dr du dx 

which may be temporally discretized to obtain a linearized representation for the source term at grid point 
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Taking the indicated derivative generates the following linear algebraic expression in terms of uf +l : 
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Substitution of this expression into the previously developed discretization equations (50) and (59) yields 
the following forms: 
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Discretization equation (52) at the symmetry node (/= 1) is unaffected by this linearization procedure and 
requires no modification. 

3.1.2 Solution Algorithm 

Now that the numerical solution has been reduced to a system of linear algebraic equations, we 
find it advantageous to adopt a linear vector space notation. Application of the preceding discretization 
equations at each nodal point, for instance, yields a system of equations that may be compactly expressed 
using the matrix equation 


KU m+l +LU m +M= 0 , (66) 

where K, L, and M are N x N square matrices containing numerical coefficients and U is an A-elcmcnt 
column vector containing the dependent variable values at the grid points. As before, superscripts denote 
temporal indices. The general solution of equation (66) at time step m + 1 follows immediately and has the 
basic form: 


U m +l =K -l LU m + M , 


(67) 


where is the matrix inverse of K. 


Equation (67) may be solved by a number of algorithms, but the most convenient and simple 
follows from the standard Gaussian elimination procedure. Because the non-zero elements of K align 
themselves along the central three diagonals of the matrix, the elimination process turns into a particularly 
simple recurrence sequence. This method is commonly referred to as the TDMA (Tri-Diagonal Matrix 
Algorithm). 30 


form 


To implement TDMA, we first proceed by writing the 
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generalized discretization equation in the 
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where 


i = 1 

i = 2,...,N — 1 (69) 
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Formal derivation of TDM A may be found in standard textbooks on numerical analysis (see 
ref. 30, for instance). The algorithm procedures are summarized as follows: 


Start the recurrence process by calculating the parameters. 


fi = - Qi 

a \ 


a l 


Continue the TDMA recurrence sequence over the range i= 2, 3, . . ., N using the relations 


(73) 
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At the other end of the sequence, note that b N = 0, which leads to P N = 0. Hence, 

m + 1 _ 

U N -QN ■ 

The remaining unknown grid values may then be determined by marching backward (/= 
N- 1, . . ., 1) using the back substitution relation 


(74) 
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3.2 Validation 

The numerical methodology for nonstationary thermal ignition was implemented in FORTRAN 
programming language to enable automated machine processing on a digital computer platfonn. As a 
validation measure, comparative baseline calculations were then performed with respect to published 
solutions in the peer-reviewed literature. For simplicity, these calculations were all performed assuming 
an initially unifonn temperature distribution within the reactive material as defined by the relation 

u(£,0) = U + C , (77) 

where U is the dimensionless ambient temperature and C is an arbitrary constant representing an initial 
perturbative displacement from thennal equilibrium. 

Given values of A, U, and Bi, our principal interest becomes a determination of the critical value 
for C, denoted by C cr , a threshold value which distinguishes between those initial conditions that converge 
to an upper steady state branch associated with ignition and those that converge to a non-igniting stable 
lower branch. Validation of the nonstationary numerical model is established through direct comparison 
of predicted values for C cr with the published results of Weber et al. 4 

Validation calculations were perfonned for all three principal centrosymmetric solids using Crank- 
Nicholson time integration (/= 0.5) with the controlling physical parameters defined as A= 10 6 , £7=0.05, 
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and Bi —> Critical initial threshold temperatures were determined through trial and error, and represen- 
tative transients for the dimensionless centerline temperature are shown in figures 13, 14, and 15 for n = 0, 
n= 1, and n= 2, respectively. Note that the centerline temperature directly corresponds to the maximum 
value in the material due to the imposed physical symmetry for the problem. 



X 

Figure 13. Centerline temperature transients for n = 0 (A = 10 6 and Bi — > °o). 


Numerically predicted threshold values for ignition were C cr =0. 0043, C cr = 0.0083, and C cr = 0.0 101 
for the planar, cylindrical, and spherical geometries, all of which are in exact agreement with the published 
results of Weber et al. 4 The stable upper steady state branches are not displayed in these figures since it is 
several orders of magnitude larger than the lower steady state branches. It should be pointed out in pass- 
ing, moreover, that the upper branch serves only as a mathematical attractor and has no practical physical 
interpretation, since ignition would occur long before the upper branch could be attained. For fundamental 
demonstration purposes, an upper steady state branch for an igniting spherical geometry case (C=0.011) 
is shown in figure 16. The ignition delay, the period of time required for the initial condition to evolve to 
the upper steady state branch, may be directly inferred from these types of graphs. 

The numerical methodology was also evaluated for finite surface convective heat transfer rates to 
validate proper response to arbitrary environmental boundary conditions. These calculations were per- 
formed for the planar slab geometry using Crank-Nicholson time integration (/= 0.5) with the controlling 
physical parameters defined as A= 10 6 and £7=0.03. The resulting variation in C cr as a function of Biot 
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Figure 14. Centerline temperature transients for n = 1 (A = 10 6 and Bi —> 



Figure 15. Centerline temperature transients for n = 2 (A = 10 6 and 5/ — > °°). 



X 

Figure 16. Upper steady state branch for n = 2 (A = 10 6 and Bi — > °°). 


number is shown in figure 17, with the exact numerical values tabulated as an inset to the graph. These 
results are also in excellent agreement with the published data of Weber et ah, 4 and demonstrate that the 
initial condition threshold for ignition becomes relatively insensitive to changes in Bi once it exceeds a 
value of about 5. 

As a final validation step, calculations were performed for a practical case involving the self-heat- 
ing of aerated milk powder, which is frequently stored in an ambient temperature environment immedi- 
ately following a high-temperature drying procedure. This particular case, based on an actual self-heating 
fire in New Zealand, 31 was also analyzed by Weber et ah, and their results are again used as a comparative 
baseline. 4 The characteristic physical parameters for milk powder are well established, and the planar slab 
geometry was considered an adequate approximation to the actual physical conditions associated with the 
incident. The established dimensionless parameters for the thennal ignition model are A=2.88 x 10 12 , 
given a half- width of 0.2 m, and A = 11.5 x 10 12 , given a half- width of 0.4 m, with U= 0.0241. Represen- 
tative transients for the dimensionless centerline temperature are shown in figures 18 and 19 yielding the 
critical threshold values C cr = 0.0058 and C cr = 0.0044 for slab depths of 0.2 m and 0.4 m, respectively. 
These results are again in excellent agreement with published results 4 

Based on the demonstrated agreement and excellent concordance with established results for a 
variety of problem conditions, it was concluded that the proposed numerical methodology had been tested 
and validated to a satisfactory extent. Although not discussed in specific detail, the methodology was also 
thoroughly exercised with both explicit and implicit source term linearization schemes, and the resulting 
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Figure 17. Variation in C cr as a function of Bi with A = 10 6 and U = 0.03. 



Figure 18. Centerline temperature transients for milk powder (A = 2.88 x 10 1 -). 
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Figure 19. Centerline temperature transients for milk powder (A = 11.5 x 10 12 ). 


numerical predictions were found to be completely insensitive to the linearization method being utilized. 
On the basis of these extensive validation efforts, it was concluded that the numerical methodology estab- 
lished in this study could serve as a reliable and dependable tool for undertaking a thorough study of the 
critical initial conditions for thennal ignition. 
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4. ASSEMBLY PROBLEM AND CRITICAL INITIAL CONDITIONS 


With the establishment of validated stationary (steady-state) and nonstationary (transient) models 
for thermal ignition/explosion, it becomes possible to undertake a judicious examination of the so-called 
‘assembly’ problem (see section 1.2), in which the fate of a self-heating material can be strongly depen- 
dent on initial conditions. Relevant physical examples include the safety of processed reactive materials 
when placed into piles or bins while hot and the initiation of explosives by localized hot spots. From a 
practical perspective, it is extremely important to know the models’ respective ranges of applicability for 
arbitrary circumstances of assembly; therefore, a deeper understanding of the underlying mathematical 
structure is highly desirable. 

The central issue at hand is whether the initial temperature excess in the assembly will continue 
to rise indefinitely or dissipate through conductive transport and Newtonian cooling to the surroundings. 
Of particular interest to this study is the critical demarcation boundary between those initial conditions 
serving as an onset to ignition/explosion and those which ultimately decay to a self-extinguishing quench 
state. For the general case of a non-uniform assembly, the transient solution for the evolved state may be 
expected to display a marked sensitivity to the spatial concentration of thermal energy at the outset, and 
criticality will depend on the initial temperature distribution. In this section, we examine such subtle math- 
ematical issues in detail by incorporating a generalized self-consistent treatment of the initial tempera- 
ture distribution for the non-stationary model, which pennits direct comparison of criticality predictions 
against the stationary model. 

The purpose is to precisely determine the conditions under which use of the steady-state model- 
ing approach is legitimate or illegitimate, and to establish deeper insight into the structure of the solution 
space. It should be noted as a warning, however, that reactants are consumed in all real systems and no true 
steady states are physically possible. As such, the temperature rise accompanying self-heating can only 
exhibit a maximum before ultimately decaying to ambient as the exothermic reaction approaches comple- 
tion and is overtaken by cooling, even after large increases associated with explosive effects, and it is no 
longer possible to recognize ignition/explosion as a fundamental discontinuity. Here, we have followed 
conventional practice and neglected reactant consumption by reasoning that if the heat of reaction is high 
enough, concentration changes preceding the onset of a large and sudden temperature rise, corresponding 
to explosion or ignition, will be negligibly small. The validity of conclusions drawn from this study is 
therefore limited by the validity of this basic underlying assumption. 

4.1 Initial Shape Profile 

As part of the previous validation exercises in section 3, we adopted the following convenient 
expression for an initially uniform temperature distribution within the material assembly: 

u(£,0) = U + C , (78) 
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where U is the dimensionless ambient temperature and C is an arbitrary constant representing the initial 
perturbative displacement from ambient conditions. The principal objective in that case, given values for 
U, A, and Bi, was determination of the critical value for C, denoted by C cr , which defines the watershed 
threshold between those initial conditions that evolve to an upper steady-state (i.e., ignition/explosion) 
and those that evolve to a self-extinguishing steady-state quench. 

We now extend this basic approach by introducing a one-parameter shape profile g(£) to define a 
generalized nonuniform initial temperature distribution having the fonn 

u(Z,0) = U + Cg(Z) . (79) 

Again, we seek critical values, C cr , such that any initial conditions with C<C cr evolve to the lower (cool) 
steady state and any conditions with C> C cr evolve to the upper (hot) steady state. The initial shape profile 
is defined by the one-parameter function: 

g(S) = Ae( 1-f). £20 

<?(<=)— * 1 5 £— >°° , (80) 

where e is the initial profile shape factor and A £ is a nonnalization factor for preserving total heat content 
irrespective of profile shape. With this definition, we anticipate the existence of a family of critical initial 
condition profiles for each value of the initial shape factor, e. 

Normalization is achieved through a constraint on the spatial integral of g(<§) such that it is invari- 
ant with e. The value of this integral is taken to be unity for convenience: 

/= Jo*(£K= i- < 81 > 

Thus, substitution of equation (80) and execution of elementary integration steps yield the following 
expression for the normalization factor: 


£+1 

A e =— , (82) 

£ 

which clearly demonstrates that A £ must decrease as £ increases in order to preserve energy content in the 
distribution. Nonnalized initial shape profiles for representative values of the shape factor are shown in 
figure 20. These include linear (£= 1), parabolic (£=2), and uniform (£=°°) initial temperature distribu- 
tions for the three principal centra symmetric solids of interest to this study. The corresponding normaliza- 
tion factors are summarized in table 4. 

4.2 Criticality in Nonuniform Assemblies 

Our principal concern is to extend the scope of analysis to nonunifonn assemblies and detennine 
the variation in C cr with U, provided there are fixed values of A and Bi, using £ as an independent param- 
eter for initial heat concentration. We may reasonably anticipate from physical intuition the following 
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Figure 20. Normalized initial shape profiles for representative values of e. 


Table 4. Summary of profile shape factors and nonnalization factors. 


Profile Shape Factor (e) 

Normalization Factor ( f E ) 

1 (linear) 

2 

2 (parabolic) 

3/2 

°° (uniform) 

1 


functional dependence: (1) C cr will decrease as U increases; (2) C cr will decrease as A increases; and (3) 
C cr will decrease as Bi decreases. It is intuitively obvious, however, that A £ C cr would better serve as the 
parameter of choice in defining the critical threshold, since A £ exhibits a strong dependence on the value of 
e. Previous cursory numerical investigation has supported the basic preceding conclusions, 4 and the intent 
here is to buttress this argument with additional computational results. An additional outcome from the 
extended numerical analysis is improved understanding and deeper insight into the mathematical structure 
of the solution sets to the model problems. 

To implement this numerical study, it was first necessary to construct an automated search proce- 
dure by which the nonstationary model could be used to determine C cr for arbitrary values of U, given 
fixed values for A, Bi, and e. This was accomplished in practice by assuming a very small value for C 
and computing a fully evolved temporal solution for a selected value of U using the coded numerical 
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methodology for the nonstationary model. If the long-term solution is found to evolve to the lower steady 
state solution, the results are recorded, C is given a small incremental increase, and the calculation is auto- 
matically repeated again and again. The value of C cr is ultimately resolved when a fully evolved solution 
to the upper steady state attractor is obtained. By starting with a very small value for U and repeating the 
above procedure with incrementally increasing values, it is possible to fully construct the critical initial 
threshold curve out to the bifurcation point, U cr . 

The above automated search procedure was implemented in computer code, and calculations were 
then carried out for the three principal centrosymmetric solids given A= 10 4 , 10 6 , 1 0 8 , assuming Bi — > °°. 
Although this crude search procedure proved effective and accurate for the purpose at hand, it is obviously 
inefficient and computationally expensive; the exact computing time being dependent, of course, on the 
size of the step increments and the speed of the microprocessor. In practice, it was necessary to balance the 
need for fidelity and resolution against the time required to construct a solution. Construction of a single 
threshold curve of satisfactory accuracy, for instance, would normally take several days of continuous run 
time on a dedicated personal computer. 

Computational results for the nonuniform assembly problem are shown in figures 21-38, which 
summarize criticality characteristics for the planar slab, infinite cylinder, and sphere using linear (e= 1), 
parabolic (e=2), and uniform (e=°°) initial temperature distributions. For each value of A and geometry 
shape, the critical threshold curves of w(0) = U+C cr g(0) versus U as well as the quench state curve «(()) 
versus U for each of the three simple geometric configurations under consideration. For purposes of com- 
parison, stable lower and unstable intennediate branch solutions of the stationary model are also displayed 
on these graphs as dashed lines. The corresponding curves for A £ C cr versus U are also displayed since they 
will prove to be of significant importance in the analysis to follow. 

The stable upper steady state curve is not shown since it is several orders of magnitude larger than 
the criticality and lower steady state branches and is of no physical consequence other than serving as a 
mathematical attractor for ignition. 

The mathematical feature of essential importance and interest to be extracted from these results 
is the fact that the exact criticality threshold, as predicted by the nonstationary model, is significantly 
different from the unstable intermediate steady state, as predicted by the stationary model. Although the 
models yield identical predictions for the bifurcation point critical ambient temperature, U cr , the ignition 
threshold for assembly of initially hot materials is found to be considerably lower than that deduced from 
the stationary model. This has serious implications since the widely held consensus view has always been 
that the unstable intermediate steady-state solution provides an adequate estimate for safe assembly of 
potentially hazardous self-heating materials. It is clear from these results, however, that this is a severely 
flawed postulate, the careless use of which could have grave practical consequences. 

A particularly interesting point of additional note is the observation that the threshold temperature 
for ignition increases as £ decreases and the initial heat distribution becomes more concentrated near the 
centerline. This seems somewhat counterintuitive, but can be explained from a heat balance perspec- 
tive. That is, the higher temperature gradient associated with more concentrated heat deposition near the 
centerline increases outward heat conduction toward the surface boundary enough to outpace the self 
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Figure 21. Criticality characteristics of nonuni form planar slab assemblies (A = 10 4 ). 


heating rate. It might be anticipated that this rate imbalance could be reversed given a sufficiently high 
exothennicity, but most practical substances of interest would not yield values very much beyond A= 10 8 . 
Thus, a uniformly distributed assembly of hot material is potentially more dangerous than a nonuniform 
assembly with a hotter interior, given the same total heat content. The shape of the assembly also has a 
notable effect on ignitability in that the critical initial threshold temperature increases as we go from planar 
to cylindrical to spherical geometry. This is partly a surface-to-volume effect and is partly caused by the 
change in heat concentration. 

Based on the preceding observations, one might speculate that the criticality threshold should 
somehow be correlated with the spatial moments of the initial heat distribution. The validity of this math- 
ematical conjecture will be taken up in section 5. 
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Figure 24. A £ C cr versus U curves for nonuniform cylindrical assemblies (A = 10 4 ). 



Figure 25. Criticality characteristics of nonuniform spherical assemblies (A = 10 4 ). 
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Figure 36. A £ C cr versus U curves for nonuniform cylindrical assemblies (A = 10 8 ). 



u 


Figure 37. Criticality characteristics of nonuniform spherical assemblies (A = 10 8 ). 






5. CORRELATION AND REDUCTION OF CRITICAL ASSEMBLY CONDITIONS 


In the previous section, an extensive numerical parametric study was undertaken to identify the 
critical threshold separating those initial conditions which lead to ignition/explosion and those which 
ultimately decay to a self-extinguishing quench state. Specifically, this study examined the influence of 
nonuniform thermal assemblies within the principal centrosymmetric solids. The outcome, a broad set of 
numerically constructed ignition thresholds over a wide range of values for the dimensionless eigenvalue, 
provided a reliable database for precisely determining ignition hazards given very specific physical con- 
ditions. From a deeper theoretical perspective, however, there is a natural desire to seek further means 
of reducing, ordering, and correlating these results such that the underlying mathematical structure may 
be revealed and expressed in the most compact and economical form possible. In this section, we pursue 
various methods of attack for achieving this objective. 

First, we demonstrate that the critical threshold curves from the parametric study can be accurately 
correlated using a general second degree conic section. We then introduce a new dimensionless parameter 
into the resulting hyperbolic correlation and develop a simple relationship that collapses the entire data set 
from the parametric study onto a single universal line. This is an entirely new result which packs accurate 
predictive power within a highly compact form. Second, we formulate a conjecture and examine spatial 
moments of the initial energy content integral as an alternative means of reducing and correlating the criti- 
cal threshold conditions. 


5.1 Correlating Forms and Structural Compaction 

At this stage, we are confronted with a large parametric database for which there are no obvious 
means of ordering and reducing the results in a straightforward deductive way. What we are seeking, of 
course, is a means of discovering and revealing hidden mathematical structure through an indirect process 
of intuitive reasoning. Under these circumstances, it is natural to examine the data in search of basic math- 
ematical fonns that can organize the results and provide improved insight. Of the numerous correlating 
forms that could be considered, conic sections offer one of the simplest and most widely useful constructs 
available. Here, we examine in detail the application of a second degree conic section correlation to the 
critical initial condition thresholds, and demonstrate effective compaction of the solution space’s math- 
ematical structure. 

5.1.1 Hyperbolic Correlation 

Close examination of the threshold curves in the {A £ C cr ,U} plane reveals structural features that 
tend to be associated with a hyperbolic conic section. To pursue the mathematical implications of this 
observation, we begin by recalling the general form of the origin centered hyperbola from elementary 
analytic geometry, as illustrated in figure 39. In this case, a point (x,y) is on the hyperbola with vertices 
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Figure 39. Origin-centered hyperbola. 
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fL_Z_ = i 
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(83) 


(±n,0) and foci (±c,0) if, and only if, it satisfies the equation where b 2 =c 2 -a 2 . Recall also that a hyper- 
bola has a pair of slant asymptotes with slope ±b/a, a distinct property that is in close correspondence to 
structural characteristics exhibited by the critical threshold curves. 


To develop the desired hyperbolic correlation for the threshold curves in the {A £ C cr ,U} plane, it 
is now necessary to translate the coordinate axes according to a well known theorem which states that a 
point (x,y) is on the hyperbola with center (h,k), vertices (h ±a,k ), and foci ( h±c,k ) if and only if it satisfies 
the equation 


(x-h) 2 (y-kf 
a 2 b 2 


(84) 


where b 2 = c 2 -a 2 . In this particular case, we identify the transverse and conjugate coordinates in the 
{A £ C cr ,U} plane as x=U and y=A £ C cr and define a new parameter r=A £ C cr to arrive at the relation 


(U~hf (r-kf 

2 ,2 

a b 


(85) 


52 



At this point, we impose the obvious constraints h=U cr +a and k= 0 and deduce our essential correlation 
equation 


(U-U cr -a) 2 r 2 

2 “72 

a b 


( 86 ) 


Utilization of this correlation equation requires the introduction of certain additional assumptions 
in order to fully define values for the vertices and the foci. As a basis for the first assumption, we note 
that the slope of the critical threshold curve r=A £ C cr away from the critical ambient temperature, U cr is 
approximately unity. Since this slope should match the slope of the slant asymptotes for the hyperbola, 
which by definition must take the value ±b/a, we immediately deduce the constraint b~—a. In this case, 
equation (86) takes the simpler form 


r 2 =(U-U cr -af -a 2 . (87) 

The final remaining issue concerns definition of the vertex, a. From casual inspection of the criti- 
cal threshold curves, we clearly anticipate that the value of a must decrease as £ increases. As such, we 
introduce the following two-parameter correlating expression: 


a =a+ m — 


( 88 ) 


where a and /3 are fitting parameters and a is assumed to vary inversely with £. Note that a multiplicative 
factor of 2 was introduced into equation (88) as a means of scaling /3 such that a = a+ f5 when £=2. 


The objective now is to deduce the best parameter values for fitting equation (87) to the precise 
critical threshold curves in the {17,(7} plane based on the complete data set from the parametric study of 
section 4. The results of this fitting exercise yielded an array of optimal fitting parameters having a one- 
to-one correspondence with the array of (n,X) values from the parametric study. These fitting parameters 
are summarized in table 5. The actual hyperbolic correlation curves are shown with the precise computa- 
tional threshold curves in figures 40-48 for each considered value of £. Figures 40-42 depict results for 
the slab geometry (n = 0) for A= 10 4 , 10 6 , 10 8 , respectively. Figures 43-45 depict results for the cylindrical 
geometry (n= 1) for A= 10 4 , 10 6 , 10 8 , respectively. Figures 46-48 depict results for the spherical geometry 
(n = 2) for A= 10 4 , 10 6 , 10 8 , respectively. 


Inspection of these figures clearly demonstrates that the hyperbolic correlation provides an excel- 
lent fit to the precise computational results over the full range of £ and A values under consideration. 
Thus, we conclude that the mathematical structure of the complete solution space is accurately captured 
by a simple hyperbolic conic section, as defined by equation (87), when using a two-parameter correlat- 
ing expression for the vertex, as given by equation (88). Variations in the fitting parameters appear to be 
smooth and well behaved, and it is clear that a and /3 could further be expressed as analytical functions of 
n and A, as well. Here, however, we shall be content to leave them in tabulated fonn, only. 
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Table 5. Optimal fitting parameters for hyperbolic correlation. 
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Figure 40. A £ C cr versus U results for nonuniform planar slab assemblies (A = 1 0 4 ). 

Precise computational value — symbols. Hyperbolic correlation — curves. 
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Figure 41. 


Figure 42. 



A £ C cr versus U results for nonuniform planar slab assemblies (A = 10 6 ). 
Precise computational value — symbols. Hyperbolic correlation — curves. 



A £ C cr versus U results for nonuniform planar slab assemblies (A = 10 8 ). 
Precise computational value — symbols. Hyperbolic correlation — curves. 
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Figure 43. 


Figure 44. 



U 


A £ C cr versus U results for nonuniform cylindrical assemblies (A = 10 4 ). 
Precise computational value — symbols. Hyperbolic correlation — curves. 



A £ C cr versus U results for nonuniform cylindrical assemblies (A = 10 6 ). 
Precise computational value — symbols. Hyperbolic correlation — curves. 
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Figure 45. A £ C cr versus U results for nonuniform cylindrical assemblies (A = 10 8 ). 

Precise computational value — symbols. Hyperbolic correlation — curves. 



Figure 46. A £ C cr versus U results for nonuniform spherical assemblies (A = 10 4 ). 

Precise computational value — symbols. Hyperbolic correlation — curves. 
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Figure 47. A £ C cr versus U results for nonuniform spherical assemblies (A = 10 6 ). 

Precise computational value — symbols. Hyperbolic correlation — curves. 



Figure 48. A £ C cr versus U results for nonunifonn spherical assemblies (A = 10 8 ). 

Precise computational value — symbols. Hyperbolic correlation — curves. 
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5.1.2 Reduction to Compact Universal Form 

Application of the hyperbolic correlation to the complete parametric solution space provided pleasing 
enough results to encourage further examination of the mathematical implications. For example, an imme- 
diate follow-on question arises as to whether it is possible to achieve additional reduction through the 
deduction of a compact universal correlating form. To pursue this question in detail, we revisit our hyper- 
bolic correlation in the simplified form of equation (87). 

Our first course of action is to define a new dimensionless temperature AU= U cr - U so that equa- 
tion (87) may be written as 


r 2 =( r AU-a) 2 -a 2 . (89) 

Then, we expand the bracketed square in order to find 

r 2 =AU 2 +2aAU (90) 

or 

AU 2 + 2aAU-T 2 =0 . (91) 

Noting that this equation is quadratic in AU, we may immediately write an expression defining AU explic- 
itly in terms of P. 


AU = \l r 2 +ci 2 -a . (92) 

We now have a convenient explicit correlating form that completely captures the essential math- 
ematical structure of the parametric solution space, and it is strikingly simple. An even simpler form can 
now be had by defining the right hand side of equation (92) as a new dimensionless criticality parameter, 

r = yjr 2 +a 2 - a . We therefore arrive at the universal linear correlating form 


r = AU . (93) 

This result is utterly simple and, in a sense, mathematically beautiful since it collapses the entire solution 
space onto a single line in the { P ,AU) plane. The results of this exercise are summarized in figures 49, 
50, and 51 for the slab, cylinder, and sphere, respectively. Although all of the data collapses to the same 
universal line, the results are depicted separately by geometry in order to avoid excess clutter. 

The fit to the universal line is generally good except for some slight deviation in the low ambi- 
ent temperature region where U«U cr This deviation arises from the fact that the actual threshold curve 
begins to inflect as it nears and passes through the U= 0 axis, whereas the hyperbolic correlation requires a 
continuous merging with the slant asymptote. Even there, the deviation is not extremely large. In general, 
the fit is better at large 8 and slightly degrades as the value becomes smaller and the initial thennal distri- 
bution has a higher degree of nonuniformity. 
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Figure 49. Universal correlating line. Slab (n = 0). 



Figure 50. Universal correlating line. Cylinder (n = 1). 
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Figure 51. Universal correlating line. Spherical (n = 2). 


5.2 Spatial Moments of Energy Content Integral 

Intuitively, we expect that the ultimate fate of a combustible assembly should have a direct math- 
ematical coupling to its initial energy distribution. We further anticipate that the predisposition for this 
assembly to either ignite or quench may be conveniently and compactly expressed as an integral over this 
initial energy distribution. Our purpose then is to rigorously define the initial energy content integral and 
precisely state a mathematical conjecture relating spatial moments of this integral to the fate of thermally 
nonuniform assemblies of combustible substances. The results of the previous computational parametric 
study may then used as a basis for thoroughly testing the conjecture and establishing a provisional basis 
for its acceptance and application to industrial fire problems of great practical interest. 

5.2.1 Definitions and Conjecture 

From physical reasoning, we assert that the critical condition for ignition of a thermally nonuni- 
form combustible assembly depends largely upon the initial energy concentration or density, particularly 
within a centralized hot region as represented by the one-parameter shape profile of figure 20. Here, the 
hot spot has been centered on the axis of symmetry as a worse case representation since the dispersion and 
removal of excess heat is most impeded under these circumstances. 


(94) 


We therefore begin our development with a fundamental definition of energy density e: 

e = pc p T , 
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where the variables on the right hand side are all as previously defined. The total energy content E must 
then follow from a volume integral over the full assembly: 


E =^>edV . 


(95) 


Now, recall our previous definitions for the nondimensional variables a=pc p /k and u=RT/E a and 
eliminate c p and T in equation (94) to obtain a dimensionless local energy density in the form: 


e = e 


akE, 

R 


a I = u 


(96) 


Pursuing a similar train of thought, we introduce the following definition for the dimensionless total 
energy content: 


E = 



and refonnulate the volume integral of equation (95) to obtain 


(97) 


E = $edV = $,udV . (98) 

Note that the total energy content within the assembly has been succinctly and conveniently expressed as 
the volume integral of the dimensionless reactant temperature over the complete assembly volume. With 
more specificity, we now apply the energy content integral to the principal centrosymmetric solids and 
reduce to the simplified expression 


1 1 

E = je(%,t)d% = ju(%,t)d% , ( 99 ) 

0 0 

which is to be evaluated at time t= 0. As a point of clarification, it should be noted that certain geometrical 
multiplication factors have been neglected here in order to achieve commonality of form. 

Let us proceed further by considering spatial moments of our dimensionless energy content inte- 
gral. The conventional approach for defining a spatial moment is to weight the kernel of the integral by the 
distance from the origin (i.e., the moment arm) raised to some arbitrary power, which defines the order of 
the moment. Naturally, this places less weight on contributions to the integral that are closer to the origin 
than those that are farther away. Because our hot spot is located on the axis of symmetry, however, we 
wish to place more weight on integral contributions that are closest to the origin. Thus, we henceforth 
define our moment arm using the factor ( 1-q) and define the /nth order spatial moment for the nth class 
geometry as 


1 

Vm,n = \{ l -£) mu {Z$) d Z • 


0 


( 100 ) 
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Recalling the generalized one-parameter initial shape profile from section 4, we deduce the practical 
working form 


1 


1 


¥m,n = J(1 “ [U + C cr g(£)] d% = J(1 - £)'” 


t/ + C cr A e (l-<f)]^ 


( 101 ) 


0 0 

It should be clear that this expression is explicitly dependent on the geometrical index n and that the influ- 
ence of geometry on the spatial moment is carried entirely within the critical constant, C cr 

We now state, based principally on intuitive reasoning, the following mathematical conjecture: 


Given any thermally nonuniform assembly in the form of a principal centrosymmetric solid, 
there exist nonstationary solutions to the Burnell-GrahamEagle-Gray-Wake reaction-diffusion 
equation, under the constraint for which certain order-/;; spatial moments of the critical 

total energy content integral, \ff n m , are invariant with changes in the dimensionless ambient tem- 
perature, U, and the initial shape profile parameter, e, within some quantifiable error band, and 
are therefore functionally dependent on the dimensionless eigenvalue, A, only. 


5.2.2 Provisional Validation 


As a means of generating reliable quantitative evidence for the previously stated conjecture, we 
now consider the first and second order spatial moments (m= 1,2) of the critical total energy content inte- 
grals over the full parameter space of the computational study. Evaluation of these integrals was effected 
by numerically integrating equation (101) using a cubic spline interpolation procedure to define C cr A e 
between available data points. The results of these evaluations are summarized in figures 52-57, which 
display the variation in y7 m , n with U/U cr over the entire parameter range of the study. 

Inspection of the resulting first- and second-order spatial moments provides provisional confirma- 
tion of our conjecture and indicates that the second-order moment, which has a narrower error band, may 
serve as a more accurate basis of prediction. It should also be noted that the error band exhibits a tendency 
to narrow as the value of A increases. Results for the third order spatial moments, which were also com- 
puted but not displayed, showed a re -widening of the error band. Apparently, the second order spatial 
moment satisfies our conjecture with the least variance and the smallest realizable error band. Spatial 
moments of noninteger order were not considered. 

If we accept these results as provisional validation for the conjecture, we now have at hand 
a simple method for rapidly calculating fire hazard risks for thermally nonuniform assemblies of self- 
heating materials. Assuming that the assembly temperature profile it (g) is known with a reasonable degree 
of confidence and that the physical characteristics of the substance are well enough defined to compute 
A to good accuracy, one may simply evaluate the spatial moment from equation (100) and compare this 
result with the mean value of the moments determined by this study. If the computed spatial moment is 
near to or greater than the mean reference value, at the specified value of A, then the risk for spontaneous 
ignition should be considered high. If the computed spatial moment is less than that mean reference value, 
one may safely assume that the risk is low or negligible, depending on the magnitude of the difference. 
Arithmetically averaged values for the spatial moments integrals i ]7 nm are summarized in table 6 for con- 
venient use and are also plotted as a function of A in figure 58. 
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Figure 52. y/ nm versus U/U cr for nonuniform slab assemblies (n = 0, m = 1). 



U/U 

cr 


Figure 53. y/ n , m versus U/U cr for nonunifonn slab assemblies (n = 0, m = 2). 
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Figure 54. y/ n m versus U/U cr for nonuniform cylindrical assemblies (n = l, m = 1). 



Figure 55. y/ n m versus U/U cr for nonuniform cylindrical assemblies (n = \, m = 2). 
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Figure 56. y/ n m versus U/U cr for nonunifonn spherical assemblies (n = 2, m = 1). 



Figure 57. y/ n m versus U/U cr for nonunifonn spherical assemblies (n = 2, m = 2). 
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Table 6. Arithmetically averaged values of spatial moment integrals, y r m n 


Slab (« = 0) 

TT 

o 

II 

A= 10 6 

A= 10 8 

m = 1 

0.041 

0.028 

0.0022 

m = 2 

0.028 

0.019 

0.015 

Cylinder (n = 1) 

rr 

O 

II 

A= 10 6 

A= 10 8 

m = 1 

0.045 

0.030 

0.022 

m = 2 

0.031 

0.021 

0.015 

Sphere (n = 2) 

o 

II 

A= 10 6 

A= 10 8 

m = 1 

0.047 

0.031 

0.023 

m = 2 

0.033 

0.021 

0.016 



Figure 58. Arithmetically averaged values of spatial moment integrals, y/ m n 
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6. DYNAMIC BOUNDARY CONDITIONS 


Frequently, self-heating assemblies of hazardous materials are inadvertently exposed to time-vary- 
ing ambient temperature conditions that can affect the critical threshold for thermal ignition/explosion. 
This might occur, for instance, if the material were stored in a location having inadequate climate control 
provisions and the diurnal temperature variations happened to be sufficiently large. It could also occur dur- 
ing the transport of a holding container if it were to be haphazardly placed near a fluctuating heat source. 
How a particular assembly would respond to such dynamic regimes of ignition is therefore a question of 
great practical importance. 

Physical situations involving dynamic boundary conditions, where the ambient temperature or 
surface heat flux has a known time dependent variation, represent a class of problems that require a for- 
mal nonstationary mathematical treatment. The standard stationary model for static regimes of ignition 
is not directly applicable to this situation, and no viable construct has ever been found that would allow 
its use for approximation purposes. In the most general case, the transient solution may be expected to 
display a marked sensitivity to ambient temperature fluctuations, and we may anticipate that this will have 
a significant impact on the threshold of criticality. The objective of this chapter is to examine the effects 
of dynamic boundary conditions on the classic ‘storage’ problem and to use the results of the nonstation- 
ary model to lay the foundations for the development of an approximate solution methodology based on 
adaptation of the standard stationary model. For clarity, this investigation will only consider the exposure 
of simple centrosymmetric geometries to periodic oscillations in ambient temperature. 

6.1 Nonstationary Model Development 

Recall that U was treated as a time dependent variable in the numerical development of the non- 
stationary model, which allows the surface boundary condition to be arbitrarily specified as a dynamic 
external driving force. Therefore, U n and LP n ' ] 1 are directly specified in equation (72, i = N) at each time 
step according to some predefined temporal variation in the ambient environment. Although this dynamic 
variation can, in principle, obey any conceivable mathematical construct, the most frequently encountered 
situation of practical physical significance involves oscillatory patterns associated with diurnal variations 
in temperature. Here, we shall confine our attention to the classic storage problem assuming an infinite 
Biot number and time-varying ambient temperature. To limit the scope of study, we restrict attention to the 
uniform initial temperature distribution only. 

6.1.1 Sinusoidal Ambient Temperature Oscillation 

The simplest construct representing periodic oscillation is the sinusoidal harmonic oscillator with 
arbitrary amplitude and period. For our purposes, this is a convenient and useful representation for investi- 
gating the influence of an oscillatory surface boundary condition, and we adopt the following simple form 
for the instantaneous ambient temperature: 
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U = U + U'sm(2nz/T + 9 ) , (102) 

where (7 is the time-averaged mean value, U' is the fluctuation amplitude, T is the period of oscillation, 
and 9 is the phase angle. Our principal interest is to examine the behavior of the critical ignition threshold 
under the influence of this dynamic boundary condition over a limited range of parameter values. 

6.1.2 Provisional Verification 

Previously, in section 3.2, we examined representative ignition threshold calculations assuming 
an initial uniform temperature distribution within the reactive material and constant ambient temperature. 
Given values for A, U, and Bi, our principal interest was to determine the critical value, C cr , and verify 
nonstationary model performance. Here, we extend these calculations to include dynamic boundary condi- 
tion effects and present provisional verification results. 

The previous validation calculations were performed for all three principal centrosymmetric sol- 
ids using Crank-Nicholson time integration (/= 0.5) with the controlling physical parameters defined as 
A=10 6 , (7 = 0.05, and Critical initial threshold temperatures were determined through trial and 

error, and representative transients for the dimensionless centerline temperature were shown in figures 13, 
14, and 15 for n=0, n= 1, and n = 2, respectively. 

In this section, we revisit calculations for the n = 0 and n = 2 cases under the action of a sinusoidal 
oscillating ambient boundary condition where a very small ambient temperature fluctuation (£/'=(). 00 1 ) is 
imposed on the mean ambient temperature ((7=0.05). The period of oscillation was T= 0.1, a value which 
is substantially less than the characteristic ignition induction time, T ind . For simplicity, we have assumed 
(p= 0 . 


These results are summarized in figures 59 and 60 for the planar and spherical geometries, respec- 
tively. Note that the critical ignition threshold values for the oscillating ambient temperature cases are less 
than those obtained using fixed ambient temperature conditions. In the planar slab, for instance, we observe 
that C cr = 0.0042 for the dynamic boundary condition compared to C cr = 0.0043 for the static boundary 
condition. Similarly for the spherical geometry, we find that C c; .=0.0100 for the dynamic boundary condi- 
tion in comparison to C cr = 0.01 10 for the static boundary condition. The significant conclusion appears to 
be that fluctuating dynamic boundary conditions generally act to reduce the ignition threshold below that 
associated with the corresponding static environment. 

This conclusion has important practical consequences for the classic storage problem when the 
environment is known to exhibit nonnegligible temporal variations in temperature. Moreover, the princi- 
pal effect of these fluctuations is to lower the ignition threshold and increase the risk for an accidental fire. 
This risk is further accentuated when the stationary model is blindly applied to the storage problem under 
environmental conditions known to have inherent dynamical variation. The purpose of the current study 
is to raise awareness of these oft neglected considerations and point the way toward the development of 
practical predictive capabilities. 
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Figure 59. Centerline temperature transients for n = 0 (A = 10 6 and Bi — > <») with sinusoidal 
oscillating ambient temperature ( U =0.05, U' = 0.001, and 7=0.1). 



Figure 60. Centerline temperature transients for n = 2 (A = 10 6 and Bi °o) with sinusoidal 
oscillating ambient temperature (C = 0.05, U = 0.001, and 7=0.1). 
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6.2 Parametric Survey 


We now extend the scope of our dynamic boundary condition survey and examine ignition charac- 
teristics over a wider range of parameter values. To avoid undue clutter at the expense of minor additional 
insight gained from completeness, we restrict attention to the planar slab geometry ( n = 0) assuming an ini- 
tial uniform temperature distribution and impose the sinusoidal ambient temperature oscillation defined by 
equation (102). Again, the controlling physical parameters are taken to be A= 1 0 6 , U =0.05, and fh— 

Before proceeding with detailed calculations of the critical ignition threshold, we first use the 
nonstationary model to compute ignition induction times, as measured from the critical initial tempera- 
ture u cr , for the fixed boundary state. The results are plotted in figure 61, which shows the variation in 
T jncj as a function of U. If the time to ignition is computed from some fixed initial temperature, we would 
expect that T ind °^ exp[\/U] (ref. 1, p. 208). However, when the time to ignition is based on the critical 
initial threshold temperature corresponding to the value of U, we find that the induction time is relatively 
small when U is much less than U cr and grows slowly until U approaches U cr , at which point it experi- 
ences explosive growth in magnitude. When U«U cr , T jnd ~().\ -0.3 , but as U — > U cr we observe that 
T /m /— >0.1 very rapidly in the neighborhood of U cr This result will appear counter-intuitive unless one 
keeps in mind the fact that the induction time defined here is based on the critical initial temperature, u cr , 
which decreases as U increases. 



Figure 61. Variation in T incj as a function of U for n = 0 (A = 10 6 and Bi — > °°). 
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Based on the preceding observations, we intuitively reason that periods of oscillation longer than 
the ignition induction time should have minimal effects on the predicted critical threshold conditions. 
Therefore, the primary range of interest for the period of oscillation extends over values less than or 
equal to T ind . In order to establish representative test cases, two specific values for the period of oscilla- 
tion were chosen for detailed investigation. These included a baseline case where 7 is on the order of the 
minimum value of T ind (i.e., 7=0.1) and an additional special case where 7 is significantly less than T ind 
(i.e., 7=0.05). We do not examine any effects associated with variable phase angle and simply take <p= 0 
to satisfy the limited purpose of this study. 

Nonstationary model calculations were then undertaken using the selected values for the period of 
oscillation with fluctuation amplitudes ranging from 0 to 0.025. The results of these survey calculations 
are summarized in figures 62 and 63, which depict the critical threshold curves w(0)= f/+C cr g(0) versus 
U as well as the stationary model bifurcation curve corresponding to the fixed ambient temperature condi- 
tion. 


There are several noteworthy points of observation that may be derived from these graphs. Most 
significantly, it is evident that the oscillating ambient temperature acts to substantially reduce the criti- 
cal ignition threshold across the entire range of mean ambient temperatures. This reduction becomes 
gradually more pronounced as the fluctuation amplitude increases. Furthermore, the degree of reduction 
is found to be highly magnified when U is small and the ignition induction time is relatively short. As 
U increases, however, the critical threshold curves are found to undergo a strong nonlinearity and are 
observed to exhibit an almost discontinuous jump if U' is sufficiently large. From a practical perspective, 
it is important to note that U cr decreases with increasing U', although the magnitude of this reduction does 
not appear to be strongly affected by changes in the period of oscillation over the parameter range of the 
study. 


6.3 Correlation 

From an applied mathematics perspective, the important questions are associated with how dynamic 
boundary conditions in the assembly problem affect the critical ambient temperature and whether any 
insight may be derived from this behavior that might provide a useful predictive capability for hazard 
analysis and risk assessment purposes. The parameter space under the action of dynamic boundary condi- 
tions is immense and essentially intractable from an analytical perspective, but some constructive progress 
may be expected from an attempt at correlation of numerical results. Let us begin by examining how the 
critical ambient temperature varies with fluctuation amplitude and period of oscillation for the special test 
case under consideration. 

Variation in U cr with If and 7 is summarized in table 7. These results are also depicted graphi- 
cally in figure 64. Note that U cr declines as the magnitude of If increases and that the reduction effect 
becomes less pronounced as the period is decreased from 7=0. 1 to 7=0.05. The greatest effect occurs for 
large values of If, as would be intuitively expected. Apparently, the longer period makes the system more 
prone to ignition by allowing the upper half cycle to exert more influence in comparison to the lower half 
cycle. That is, criticality seems to be affected more by prolonged heat loss suppression during the upper 
half cycle than by prolonged heat loss enhancement during the lower half cycle. This result should not 
be blindly accepted as a definitive deduction, however, since it is impossible to draw general conclusions 
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Figure 62. Criticality characteristics (n = 0, A= 1 0 6 , U =0.05, and with sinusoidal 

oscillating ambient temperature (7=0.1). 



Figure 63. Criticality characteristics ( n = 0, A= 1 0 6 , U= 0.05, and Bi —> «) with sinusoidal 
oscillating ambient temperature (7=0.1). 
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Table 7. Variation in U cr with U and T. 


If 

© 

II 

St. 

Si* 

Z7 cr (r = 0.05) 

0 

0.05033 

0.05033 

0.001 

0.05028 

0.05028 

0.005 

0.05022 

0.05022 

0.010 

0.05010 

0.05019 

0.015 

0.04991 

0.05006 

0.020 

0.04951 

0.04991 

0.025 

0.04903 

0.04971 



U' 


Figure 64. Variation in U cr with U and T lor n = 0 (X = 10 6 and Bi — > ») with sinusoidal 
oscillating ambient temperature. 


from this preliminary and cursory analysis. Ultimately, more in-depth analysis and study will be needed to 
fully resolve such issues to reach reliable deductions. 

The question that needs to be addressed is how such results can be used for predictive purposes. Is 
there some way, for instance, that the stationary model can be modified and effectively utilized to generate 
useful predictions? Here, we suggest one simple approach based on the introduction of an effective Biot 
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number to represent the effect of oscillating ambient temperature. To do so, Bi is varied in the stationary 
model in order to obtain a value for U cr which is in agreement with the tabulated values of table 7 and 
thereby deduce an effective Biot number, Bi e jj. The results of this correlating exercise are summarized in 
the graph of figure 65, which depicts Bi e jj as a function of If for both oscillation periods of interest. The 
accuracy of this correlation is quite good and opens the door to a potentially simple predictive methodol- 
ogy for the dynamic boundary storage problem. A critical assessment of feasibility and comprehensive 
development of the approach is beyond the current scope of study and is left for future examination. 



Figure 65. Stationary model correlation for Bi e ^ as a function of If and T for n = 0 
(A = 10 6 and Bi — > oo ) with sinusoidal oscillating ambient temperature. 
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7. CONCLUDING REMARKS 


Prior to this work, a number of practically important questions remained unanswered concerning 
the critical threshold for thermal ignition, given an assembly problem with a nonunifonn initial tempera- 
ture distribution. To seek answers to such questions, the reaction-diffusion equation in the dimensionless 
form d t u=V 2 u+Xe~ l/u over the bounded region Q of the principal centrosymmetric solids was critically 
reexamined from the perspective of modern numerical methodologies. The purpose was to broadly and 
deeply penetrate the solution space in order to reveal a more detailed description of the underlying math- 
ematical structure, which could be further used to establish conjectural and correlating principles of gen- 
eral predictive utility. 

As part of this undertaking, the classic stationary model was first revisited, and an innovative solu- 
tion methodology was developed whereby the two-point boundary value problem could be reexpressed 
as an initial value problem for a system of first-order ordinary differential equations. Careful validation 
calculations demonstrate that this simplified approach was capable of yielding bifurcation diagrams just 
as accurate as more sophisticated path following techniques. 

A numerical procedure was then implemented for solving the full time-dependent form of the 
reaction diffusion equation and establishing a nonstationary thennal ignition model. This numerical devel- 
opment combined second-order central differencing for the spatial derivative with a generalized time 
integration scheme, including a time-dependent convective boundary condition, to produce an efficient 
and accurate solver routine. The resulting methodology was successfully validated against published 
solutions in the peer-reviewed literature and found robust and reliable. 

Having developed and validated stationary and nonstationary models, a judicious examination of 
the assembly problem was undertaken to expose how the fate of a self-heating material is dependent on 
the internal spatial concentration of thennal energy at the time of assembly. By introducing a nonnalized 
shape profile for the initial temperature distribution, as defined by a single geometric parameter, it was 
then possible to compute criticality threshold characteristics over a broad range of practical values for the 
dimensionless eigenvalue parameter. 

It was shown how the resulting mathematical structure for the ignition threshold curves could be 
correlated by a hyperbolic conic section with a high degree of fidelity over the full range of positive ambi- 
ent temperature values. Moreover, the clever introduction of new dimensionless parameters within this 
hyperbolic correlating form was found to generate further simplification leading to a universal correlating 
form capable of collapsing the entire solution space onto a single line in the plane of the new variables. 
This result was found to hold over a wide range of shape parameters and is therefore believed to be of 
general significance. 

One of the major open questions addressed by this study of the assembly problem concerned the 
widely held conviction that certain spatial moments of the initial temperature profile ought to possess a 
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direct mathematical link to the critical ignition threshold. As such, the /nth order spatial moment of the 
critical total energy content integrals was defined, and a conjecture was formulated stating that certain 
orders of this moment would remain invariant with changes in ambient temperature and initial shape 
profile and would therefore be functionally dependent on the dimensionless eigenvalue only, within some 
quantifiable error band. Evaluations of the spatial moment integrals demonstrated this conjecture to be 
provisionally valid, with the best accuracy obtained for second order moments. The invariance property 
for the spatial moment turns out to be quite powerful and yields a simple but fairly accurate method for 
estimating fire hazard risks for thermally nonunifonn assemblies of self-heating materials. The framing 
and verification of this conjecture represents a significant accomplishment and greatly expands our math- 
ematical understanding of this practically important problem. 

As a final focal point for this work, an attempt was made to undertake a preliminary and cursory 
analysis of the classic storage problem with dynamic boundary conditions assuming infinite Biot num- 
ber. Specifically, an examination was made of ignition characteristic for the planar slab geometry with 
initially unifonn temperature and sinusoidal oscillating ambient temperature. The critical threshold limits 
were computed assuming a range of fluctuation amplitudes for two representative periods of oscillation. 
It was found that an oscillating boundary condition generally acts to reduce the critical ignition threshold 
with the degree of reduction increasing with increasing oscillation amplitude. It was also shown how the 
stationary model could be used to correlate variations in the critical ambient temperature with changes in 
the fluctuation amplitude by defining an effective Biot number to represent a virtual heat loss suppression 
mechanism. That is, there is an equivalent finite effective Biot number for which the stationary model 
accurately predicts the critical ambient temperature corresponding to the dynamic boundary storage prob- 
lem with infinite Biot number. It should be emphasized that this cursory result, based on a limited explora- 
tion of an immense solution space, must be viewed as preliminary in nature and not a statement of general 
validity. However, it is believed that this result is highly suggestive of a potentially promising approach to 
an otherwise intractable problem and is deserving of further study and investigation. 
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